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We employ a first-principles lattice-gas Hamiltonian (LGH) approach to determine the lateral 
interactions between O atoms adsorbed on the Pd(100) surface. With these interactions we ob- 
tain an ordering behavior at low coverage that is in quantitative agreement with experimental 
data. Uncertainties in the approach arise from the finite LGH expansion and from the approximate 
exchange-correlation (xc) functional underlying the employed density-functional theory energetics. 
We carefully scrutinize these uncertainties and conclude that they primarily affect the on-site energy, 
which rationalizes the agreement with the experimental critical temperatures for the order-disorder 
transition. We also investigate the validity of the frequently applied assumption that the ordering 
energies can be represented by a sum of pair terms. Restricting our LGH expansion to just pairwise 
lateral interactions, we find that this results in effective interactions which contain spurious contri- 
butions that are of equal size, if not larger than any of the uncertainties e.g. due to the approximate 
xc functional. 

PACS numbers: 68.43.Bc, 68.43.De 



I. INTRODUCTION 

Lateral interactions between species adsorbed at solid 
surfaces are crucial microscopic quantities that have been 
the target of surface science studies for a long time. 1 
These interactions govern both the equilibrium, as well 
as the non-equilibrium ordering behavior of the adsor- 
bates, and thereby critically influence the surface func- 
tion and properties in important applications like het- 
erogeneous catalysis. Traditionally, considerable efforts 
have been devoted to determine lateral interactions em- 
pirically from experimental data, e.g. from temperature 
programmed desorption or low energy electron diffrac- 
tion (LEED) measurements. In order to simplify the in- 
herently indirect determination from sparse experimental 
data, the assumption of exclusively pairwise interactions 
between the adsorbed species has often been applied. In 
recent years, algorithmic advances and increased compu- 
tational power have made it possible to determine the lat- 
eral interactions alternatively from first-principles. Most 
notably, these are approaches that parameterize lattice- 
gas Hamiltonians (LGHs) with density-functional theory 
(DFT) energetics, also called cluster expansions .2iM^ 
Since the accuracy of the determined lateral interactions 
should be of the order of UbT to properly describe the 
thermal ordering, a concern with this approach has been 
whether the employed first-principles energetics is actu- 
ally accurate enough. 

Within this context, the present work has a method- 
ological and a materials science motivation. The method- 
ological motivation is to scrutinize both the assumption 
of exclusively pairwise interactions, and the accuracy 
with which the first-principles LGH approach can provide 
the lateral interactions. For this purpose we concentrate 
on a simple model case, namely the on-surface ordering of 
atomic adsorbates at a (100) cubic surface, for which ex- 
tensive studies with model interactions have already been 



performed ,g i L§i2 a 10 i 11 i 12 i 13 To make contact with a spe- 
cific material and with experiment, we specifically choose 
the on-surface adsorption of oxygen at Pd(100), for which 
detailed experimental data on the ordering behavior is 
available^. Since in this system higher oxygen coverages 
above ~ 0.5 monolayers [ML, defined with respect to 
the number of Pd atoms in one layer of Pd(100)] induce 
structures containing incorporated oxyge n 14 i 15 ' 16 i 17 ' 18 , 
we concentrate on the low coverage regime. For this 
regime, two ordered structures have hitherto been charac- 
terized experimentallyi^±&±i±&2°.: A p(2 x 2)-0 struc- 
ture at 0.25 ML and a c(2 x 2)-0 structure at 0.5 ML, 
both with O adsorbed in the on-surface fourfold hol- 
low sites. The material science motivation of our first- 
principles LGH study is then to extract the lateral inter- 
actions operating between the adsorbed O atoms at the 
surface and to study the ordering behavior they imply. 
Specifically, this is to see whether we can confirm the 
experimentally determined ordered structures, as well as 
the critical temperatures for the order-disorder transition 
in the low coverage regime. 

Presenting a systematic first-principles lattice-gas 
Hamiltonian expansion, we indeed find the calculated set 
of lateral interaction energies to be fully consistent with 
the experimentally reported low coverage phase diagram. 
Critically discussing the uncertainties of our approach, 
both with respect to the employed LGH expansion and 
the underlying DFT energetics, we conclude that they 
primarily affect the on-site energy. The lateral interac- 
tion energies, on the other hand, can be determined with 
quite high accuracy, which we estimate for the present 
system to be around 60meV. Comparing these interac- 
tion energies with those determined previously empiri- 
cally and using the pairwise interaction approximation 
demonstrates that the latter assumption introduces an 
error that is at least as large as this remaining uncer- 
tainty when carefully determining the lateral interactions 
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from present-day first-principles calculations. 



II. THEORY 

A. Lattice-Gas Hamiltonian 

In order to describe the site-specific adsorption of oxy- 
gen atoms on the Pd(100) surface we employ the concept 
of a lattice-gas Hamiltonian, in which any system state 
is defined by the occupation of sites in a lattice, and 
the total free binding energy of any configuration is ex- 
panded into a sum of discrete interactions between the 
lattice sites (see e.g. Refs. 0HI3I1)- For a one component 
system with only one site type, this energy reads (with 
obvious generalizations to multi-site cases) 

r 

NF b = F°"- Sitc ]>> + E^ P E n ^ + 

i u=l {i<j) u 

q 
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where the site occupation numbers ni = or 1 indicate 
whether site I in the lattice is empty or occupied, with 
a total of N sites occupied, and F ™ _sltc j s the free en- 
ergy of an isolated species at the lattice site, including 
static and vibrational contributions. There are r pair 
interactions with two-body (or pair) interaction energies 
V£ between species at nth nearest neighbor sites, and q 
trio interactions with V* three-body interaction energies. 
The sum labels (i < j) u [and (i < j < k) u ] indicate that 
the sums run over all pairs of sites (ij) (and three sites 
(ijk)) that are separated by u lattice constants, and the 
summation is done such, that each pair (or trio) of occu- 
pied sites contributes exactly once to the lattice energyiSi 
Formally, higher and higher order interaction terms 
(quattro, quinto,...) follow in the infinite expansion of eq. 
(fT|). In practice, the expansion must (and can) be trun- 
cated after a finite number of terms. Obviously, the judi- 
cious choice of which interactions to consider and which 
ones to neglect must critically affect the reliability of the 
entire approach. To quantify the impact of this choice 
on accuracy, we rely on the concept of leave-one-out cross 
validation (LOO-CV) detailed below to identify the most 
important interactions out of a larger pool of possible 
interactions. Figure [1] illustrates the lateral interactions 
contained in this pool, which range from pair interactions 
up to the fifth nearest neighbor, via all trio interactions 
up to second nearest neighbor, to several compact quat- 
tro and one quinto interaction. The pool focuses thus 
on short- to medium-ranged interactions. Interactions at 
larger distances can be substrate-mediated elastic or of 
electronic origin^, but for the present system we do not 
expect such interactions to play a role on the accuracy 
level of interest to this study. 




FIG. 1: (Color online) Top view of the Pd(100) surface, illus- 
trating the considered pool of 17 lateral interactions between 
O atoms in on-surface hollow sites. (m — 1,2,3,4,5) 

are the two-body (or pair) interactions at first, second, third, 
fourth and fifth nearest neighbor distance. V^, V£, V„ u are 
considered compact trio, quattro and quinto interactions, re- 
spectively. Light grey (yellow) spheres represent Pd atoms, 
and small dark grey (red) spheres O atoms. 



B. Static and vibrational average binding energy 

In order to generate a quantitatively accurate LGH, 
we parameterize the unknown lateral interaction ener- 
gies contained in the LGH by first-principles calculations. 
The central quantities required for this parameterization 
are computed average free binding energies for a set of 
ordered configurations of O adsorbed at Pd(100). We 
write this average free binding energy as 



F b (T) = E b + F b vih -(T) 



(2) 



separating the total and vibrational contributions, E b 
and F^ lh -(T) respectively. The former is defined as 
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gies of the surface containing oxygen, of the correspond- 
ing clean Pd(100) surface, and of an isolated oxygen 
molecule, respectively. Since a free O2 molecule is thus 
used as the zero reference for Eb, a positive binding en- 
ergy indicates that the dissociative adsorption of O2 is 
exothermic at T = OK. 

In order to determine the vibrational contribution to 
the average free binding we use the phonon density of 
states u{ui) and write the vibrational free energy as 



F vib -(T) = J duF{T,u) o{lu) 



(4) 
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is the vibrational free energy of an harmonic oscillator 
of frequency u>. 23 ks is the Boltzmann constant, and 
[3 = l/(fcgT) the inverse temperature. The vibrational 
contribution to the average binding energy can then be 
written in exactly the same way as eq. (|3|), namely 
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To evaluate this contribution in practice one must thus 
determine the difference of the surface phonon density 
of states of the adsorbate covered and of the clean sur- 
face, 0o/pd(ioo)( w ) anc ^ <J Pd(ioo) ( w ) respectively, as well 
as the vibrational frequencies of the gas phase molecule 
contained in oo 2 (gas) ( w )- 



C. Total Energy Calculations 

The total energies required to evaluate eq. ([3]) are 
obtained by DFT calculations within the highly accu- 
rate full-potential augmented plane wave + local orbitals 
(LAPW/APW+lo) scheme^, using the generalized gra- 
dient approximation (GGA-PBE) 25 for the exchange- 
correlation functional. All surface structures are modeled 
in a supercell geometry, employing fully-relaxed symmet- 
ric slabs (with O adsorption on both sides of the slab) 
consisting of five (100) Pd layers with an optimized bulk 
lattice constant of a = 3.947 A (neglecting bulk zero- 
point vibrations). A vacuum region of > 10 A ensures the 
decoupling of consecutive slabs. The LAPW/APW+lo 
basis set parameters arc listed as follows: Muffin tin 
spheres for Pd and O are i?MT = 2.1 bohr and R^t = 
1.1 bohr, respectively, the wave function expansion inside 
the muffin tin spheres is done up to Z™^ x = 12, and the 
potential expansion up to = 6- The energy cutoff 
for the plane wave representation in the interstitial re- 
gion between the muffin tin spheres is E^ l ax = 20 Ry for 
the wave functions and = 196 Ry for the potential. 

Monkhorst-Pack (MP) grids are used for the Brillouin 
zone integrations. Specifically, we use a (12 x 12 x 1) grid 
for the calculation of (1 x 1) surface unit-cells. For the 
larger surface cells, care is taken to keep the reciprocal 
space point sampling identical by appropriately reducing 
the employed k-meshes. 



To obtain the total energy of the isolated O2 molecule, 
s exploit the relation Eg*^ = 2E%$ om) - D, where 

total 



i?o(atom) is the total energy of an isolated oxygen atom, 
and D the theoretical O2 binding energy. The isolated O 
atom is then calculated spin-polarized, inside a rectangu- 
lar cell of side lengths (12 x 13 x 14) bohr, T-point sam- 
pling of the Brillouin zone and without spherically aver- 
aging the electron density in the open valence shell. For 
D we employ the previously published ultra-converged 
GGA-PBE value of 6.21 eV.^ 

For the calculations of the adsorbate vibrational 
modes, the dynamical matrix is set up by displacing the 
O atom from its equilibrium position in 0.05 A steps. An- 
ticipating a good decoupling of the vibrational modes due 
to the large mass difference between Pd and O, the posi- 
tions of all atoms in the substrate below the adsorption 
site are kept fixed in these calculations. The frequen- 
cies and normal modes are then obtained by subsequent 
diagonalization of the dynamic matrix. 



D. Monte Carlo simulations 

Once a reliable set of interactions has been established, 
evaluating the LGH for any configuration on the lattice 
corresponds merely to performing an algebraic sum over 
a finite number of terms, cf. eq. ([T]). Due to this sim- 
plicity, the LGH can be employed to evaluate the sys- 
tem partition function. Here this is done by canonical 
Monte Carlo (MC) simulations for O coverages up to 
= 0.5 ML. The employed lattice size was (40 x 40) with 
periodic boundary conditions. Metropolis sampling used 
2000 MC passes per lattice site for equilibration, followed 
by 10000 MC passes per site for averaging the thermo- 
dynamic functions. Increasing any of these numerical 
parameters led to identical results on the accuracy level 
of interest to this study, i.e. here primarily critical tem- 
peratures that are converged to within 5-10 K. 

For fixed coverage on the surface, ordered structures 
are identified by evaluating order parameters sensitive to 
lateral periodicities. To check on the (2 x 2) periodicity of 
the two experimentally characterized ordered structures, 
we divide the (100) cubic lattice into four interpenetrat- 
ing sub-lattices a, b, c and d in a (2 x 2) arrangement 
This allows to separately evaluate in the MC runs 
the average number of occupied sites in each sub-lattice, 
N a , Nt, N c , and N d , respectively. Using Fourier theory, 
the order parameter for the p{2 x 2) structure is then 
defined as 
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where Aq - N a +N b + N c -N d , N 2 = N a + N b -N c + N d , 
N 3 = N a ~ N b + N c + N d , and N tot is the total number 
of sites in the simulation cell. In the same way, the order 
parameter for the c(2 x 2) structure is defined as 
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TABLE I: Calculated binding energies Eb (in meV/O atom) 
for O adsorption in on-surface hollow or bridge sites. The 
nomenclature and geometric arrangement in the surface unit- 
cell for the five ordered adlayers are explained in Fig. [2] 

(3 x 3)-0 p(2 x 2)-0 c(2 x 2)-0 (2 x 2)-30 (1 x l)-0 

9 0.11 ML 0.25 ML 0.50 ML 0.75 ML 1.00 ML 

hollow 1249 1348 1069 643 344 

bridge 1024 961 801 573 378 



Computing these order parameters as a function of tem- 
perature, the critical temperature for the order-disorder 
transition is defined by the inflection point where 4*^(2 x 2) 
or $ c (2x2) go to zero. In parallel, we also derive the crit- 
ical temperature from evaluating the specific heat, ob- 
taining values that are identical to within 10 K with those 
inferred from the order parameters. 



III. FIRST-PRINCIPLES LATTICE-GAS 
HAMILTONIAN FOR O AT PD(IOO) 

A. Energetics of on-surface adsorption 

Owing to the tendency of oxygen atoms to prefer highly 
coordinated binding sites at late transition metal sur- 
faces, the high-symmetry fourfold hollow sites appear as 
the most likely adsorption sites at Pd(100). On the other 
hand, one cannot exclude a priori that the twofold bridge 
sites are not also metastable, i.e. local minima of the po- 
tential energy surface. To test this, we slightly displaced 
a bridge site O adatom in ap(2x 1) configuration laterally 
towards a neighboring hollow site. The resulting forces 
relaxed the adatom back to the ideal bridge position, so 
that at least in this configuration the bridge site is not 
just a mere transition state, i.e. a saddle point of the 
potential energy surface. As this might also be true for 
bridge site adsorption in other (local) O adatom arrange- 
ments, we calculated the binding energetics of O atoms 
in the fourfold hollow and in the twofold bridge sites for 5 
different ordered overlayers spanning the coverage range 
up to one ML. The periodicities of these overlayers are 
explained for the case of hollow site adsorption in Fig. [2l 
and Table |T] summarizes the calculated binding energies. 

For lower coverages the fourfold hollow site is energet- 
ically clearly more stable, which suggests an insignificant 
contribution of bridge sites to the ordering behavior at 
coverages up to around 0.5 ML, even if the latter are al- 
ways metastable sites. Although the reversal of the ener- 
getic order between hollow and bridge sites at O = 1 ML 
seen in Table U is intriguing, it clearly occurs in a cover- 
age range where surface oxide formation and eventually 
three-dimensional oxide cluster growth takes place. 17 ! 18 
Since our interest lies in the on-surface ordering behavior 
at low coverages, we will therefore focus the LGH expan- 
sion for the moment exclusively on adsorption into the 
fourfold hollow sites, and return to the role of bridge site 



O atoms in Section HID. 

As the next step in the LGH parameterization, aver- 
age binding energies for different ordered configurations 
with O atoms in on-surface hollow sites and with surface 
unit-cells up to (3 x 3) were correspondingly computed. 
Despite our focus on the lower coverage regime, the set 
does comprise structures with coverages up to 1 ML, since 
these structures are required during the LGH parameteri- 
zation to determine in particular the higher-order many- 
body interactions occurring in (locally) denser adatom 
arrangements. In most cases configurations that we ini- 
tially prepared with high coverage of on-surface O atoms 
were truly metastable, in the sense that they relaxed into 
geometries where all O atoms remained in the on-surface 
hollow sites. However, some structures directly relaxed 
into geometries with O incorporated below the first Pd 
layer. Since the 0-0 interaction in these structures does 
not correspond to the physical situation we want to de- 
scribe (and would thus mess up the chosen LGH), we ex- 
cluded these configurations from our set (using a strongly 
enlarged Pd first interlayer distance as criterion). This 
left a set of 25 configurations with on-surface O atoms, 
which was subsequently used in the parameterization of 
the LGH. The binding energy data of this set is com- 
piled in Fig. [3l and indicates already overall strongly 
repulsive lateral interactions, which reduce the binding 
energy with increasing coverage by up to 1 eV. 



B. Lateral interactions 

In order to determine the lateral interaction energies 
in the LGH, we employ eq. (Q]) to write down the LGH 
expression for each of the ordered configurations calcu- 
lated by DFT (including the interactions with the pe- 
riodic images in the neighboring cells). Neglecting the 
vibrational contributions in eq. @ for the moment, we 
equate the right hand side of eq. ([TJ) with NoE^ for 
the corresponding configurations and arrive at a system 
of linear equations that can be solved for the unknown 
values of the lateral interaction energies. The crucial as- 
pects of this procedure are therefore i) the number and 
type of interactions included in the LGH expansion, and 
ii) the number and type of ordered configurations that 
are computed with DFT to determine the values of these 
interaction energies. In the following we show how i) 
is addressed by leave-one-out cross-validation, and ii) is 
aided by a search for the LGH "ground state" structures 
and an iterative refinement of the input structure set. 

1. Leave-one-out cross-validation (LOO-CV) 

In a truncated LGH expansion with finite ranged in- 
teractions, sparse configurations will exhibit a lattice en- 
ergy that is simply iVo times the on-site energy, as soon 
as all adsorbed species have distances from each other 
that exceed the longest range considered interaction. We 



FIG. 2: (Color online) Top view of 5 ordered adlayers with in on-surface hollow sites. The coverage of each configuration 
from left to right panel is 1/9, 1/4, 1/2, 3/4 and 1 ML, respectively. Light grey (yellow) spheres represent Pd atoms, small 
dark grey (red) spheres O atoms, and the black lines indicate the surface unit-cells. 
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FIG. 3: (Color online) Coverage (O) dependence of the calcu- 
lated DFT binding energies of 25 ordered configurations with 
O atoms in on-surface hollow sites. 



therefore fix the on-site energy £:° n_slte i n e q. ([T]) to be 
the DFT binding energy computed for 1/9 ML coverage, 
cf. Fig. [2]and Table Q] In this particular (3 x 3)-0 con- 
figuration, the minimum distance between O adatoms is 
8.37 A, i.e. six nearest neighbor sites away. This is larger 
than the farthest reaching interaction contained in our 
pool of lateral interactions, cf. Fig. [TJ so that fixing 
^on-site tQ thc i/ 9 ML ( 3 x 3 )_q binding energy should 

prevent fitting noise into this parameter. 

To get some guidance as to which could be the leading 
lateral interactions to be included in the LGH expansion, 
we estimate the predictive power of the LGH by the con- 
cept of leave-one-out cross- validatio n 4 i 27 i 28 ' 29 i 30 . For a 
given set of LGH interactions the cross-validation (CV) 
score is calculated as 



CV = 



\ 



1 M 
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(9) 



Here, the sum runs over the remaining M = 24 ordered 
configurations that were calculated with DFT (apart 
from the 1/9 ML coverage structure used already for 
E° n ~ sitc ), and which have DFT calculated binding en- 
ergies E® FT (i). The quantity GH (i) of the ith config- 
uration, on the other hand, is evaluated from the LGH 
expression for this configuration, cf. eq. fT]). where the 



values of the considered lateral interactions are obtained 
from least-squares fitting to the DFT energies of the re- 
maining M — 1 calculated configurations, i.e. leaving 
exactly the ith configuration out of the fit. This way, 
the CV score is intended to be a measure of the predic- 
tive power of a LGH expansion considering a given set 
of lateral interactions. In general, one would expect sets 
containing too few interactions to be too inflexible and 
thus leading to a high CV score, whereas sets containing 
too many interactions as loosing their predictive power 
through overfitting and thereby also leading to a high 
CV score. The hope is thus to identify the optimum set 
of considered interactions as that set that minimizes the 
CV score. 

Within this approach, we evaluate the CV score for 
any set of interactions out of the larger pool of 17 lat- 
eral interactions shown in Fig. [T] Table |TT] summarizes 
these scores subdivided into the optimum sets containing 
m lateral interactions, i.e. listed are the sets that yield 
the lowest CV score for any arbitrary combination of m 
lateral interactions out of the total pool of 17. For these 
sets, we then determine the values of the considered lat- 
eral interactions by least-squares fitting to the computed 
DFT binding energies of all M = 24 ordered configura- 
tions and also include them in Table [TTJ The minimum 
CV score reached indeed decreases initially upon adding 
more interactions to the set, and then increases again 
for sets containing more than 10 interactions. Another 
gratifying feature is that almost always the same inter- 
actions are picked out of the pool, i.e. the optimum set 
for m + 1 interactions corresponds mostly to those in- 
teractions already contained in the optimum set for m 
interactions plus one additional one. Only very rarely is 
an interaction that is contained in the optimum m set 
not selected in the optimum m + 1 set. And if this hap- 
pens, this concerns lateral interactions for which only 
very small values are determined, and which are thus 
anyway not meaningful within the uncertainties of our 
approach. Also in a physical picture, the determined val- 
ues for the lateral interactions appear quite plausible for 
m up to 11. The pairwise interactions decrease with in- 
creasing distance, and the leading higher-order trio and 
quattro interactions are smaller in size than the most 
dominant nearest-neighbor pair interaction. The quinto 
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TABLE II: List of the sets containing m lateral interactions, together with their CV scores and the values determined for the 
lateral interaction energies (no entry at a position in the table means that this interaction is not contained in the set. Lateral 
interactions shown in Fig. [T] but not shown here are never selected out of the pool). Negative values for the interaction 
energies indicate repulsion, positive values attraction. The sets shown are those that minimize the CV score among all possible 
sets containing m lateral interactions out of the pool of 17 shown in Fig. [T] The first line for each set corresponds to the 
data obtained by fitting to 24 ordered DFT configurations, while the second line is obtained after adding 2 additional DFT 
configurations to the fit (see text). Units for the CV score and lateral interactions are meV. 
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interaction contained in the pool of possible interactions 
is never selected. In contrast, the m = 12 set shows al- 
ready clear signs of linear dependencies, with some trio 
interactions suddenly exhibiting very large values. This 
continues for expansions containing even more interac- 
tions (not listed in Table|T|, which exhibit more and more 
obviously meaningless lateral interaction energies. 



Another equally important feature of the expansions 
up to m = 11 is the stability of the determined lateral 
interaction values against adding further interactions to 
the set. In particular for the optimum sets close to the 
overall CV minimum, i.e. for m equal to 9 or 10, adding 
another lateral interaction to the set changes the values 
for the dominant interactions by less than 2 meV. A sim- 
ilar behavior is obtained for another test to which we 
subject our expansion: Having calculated another two 
DFT configurations, we increased the set of DFT con- 
figurations employed in the fit to M — 26 and repeated 
the entire CV score evaluation. The results are also in- 
cluded in Table |TT] and show only minimal variations for 
the sets up to m = 11. Almost always the same lateral 
interactions are picked out of the pool and also their val- 
ues change by less than lOmeV compared to the previous 
procedure employing 24 DFT configurations. For the set 
to = 12, adding the two DFT input energies removes 
the linear dependencies and brings the set back in line 
with the other sets with fewer interactions. These find- 
ings suggest that the expansions are also stable against 
adding further DFT configurations and we finally iden- 
tify the set containing 9 lateral interactions and using 
24 DFT configurations to determine their values as our 
optimum LGH expansion. 



2. Ground state search 

Before moving to the ordering behavior at finite tem- 
peratures, a crucial test for the validity of the LGH ex- 
pansion is that it gives the correct ordered ground states 
at T — K, i.e. the lowest-energy structures at a given 
coverage. Here, this refers in particular to the ground 
states predicted by the DFT energetics, since the lat- 
ter is the input with which the LGH expansion must be 
consistent. Obviously, if the energetic order of compet- 
ing configurations is wrong at the DFT level (e.g. due 
to the employed approximate xc-functional) , there is no 
hope that a correct LGH expansion could cure this prob- 
lem. To this end, it is useful to rcplot the DFT database 
shown in Fig. [3]in form of a formation energy ploti^. For- 
mation energies AEf are in general defined as an excess 
energy with respect to the equivalent amounts of pure 
constituents. For the present case of on-surface O ad- 
sorption in the coverage range below 1 ML we therefore 
define 



AE f 



1 



N, 



P total 
^O/Pd(100) 



(1 - 0)25** 



total 
(100)" 



r\ 771 total 

^'- C/ (lxl)-O/Pd(100) 



(10) 



As in eq. ([3]), ^o/Pd(ioo) is the total energy for a spe- 
cific adsorbate configuration with A^o O atoms per sur- 
face unit-cell (corresponding to a coverage = No/N to t 
with iVtot the number of sites per surface unit-cell), 
2?Pd(ioo) is the total energy of the clean surface, and 
^(ixi)-c>/Pd(ioo) ^ s total energy of the full monolayer 
(1 x l)-0 configuration. With this definition, AEf re- 
flects the relative stability of a particular configuration 
with respect to phase separation into a fraction of the 
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© (ML) 

FIG. 4: (Color online) Formation energies AEf as computed 
with DFT for the 25 ordered configurations shown in Fig. [5] 
(big (red) circles), as well as for all configurations obtained 
by direct enumeration and using the LGH expansion (small 
(black) circles), see text. The red line is the convex hull for 
the DFT data, identifying three ordered ground states shown 
as insets: p(2 x 2)-0 at 9 = 0.25 ML (left inset), c(2 x 2)-0 
at 6 = 0.5 ML (top inset), and (2 x 2)-30 at 9 = 0.75 ML 
(right inset). Pd = large, light (yellow) spheres, O = small, 
dark (red) spheres. 

full monolayer configuration and a fraction (1 — 0) of 
clean surface, and we can relate it to the binding energy 
of the configuration by 

AEf = [-E&,o/Pd(100) _ -S&,(lxl)-O/Pd(100)] ■ (11) 

Plotting AEt versus coverage as done in Fig. fallows 
to identify the ground states, i.e. lowest-energy ordered 
phases, that are predicted by the present DFT data set, 
as those lying on the convex hull (or so-called ground 
state line)^ Any structure that yields a formation en- 
ergy that is higher than this ground state line is unstable 
against decomposition into the two ordered configura- 
tions represented by the two closest lying corner points 
on the convex hull. As apparent from Fig. U the convex 
hull formed by the DFT data set exhibits three ordered 
ground states (apart from the trivial ones at the ends of 
the considered coverage range). Consistent with exist- 
ing experimental data, these are the p(2 x 2)-0 ordered 
phase at = 0.25 ML, and the c(2 x 2)-0 ordered phase 
at = 0.5 ML. A third ordered structure, (2 x 2)-30 at 
= 0.75 ML, is at best metastable, since it falls already 
in the coverage range above ~ 0.5 ML, for which surface 
oxide formation sets in. 

Using eq. (|11[) we can also evaluate formation ener- 
gies using the binding energies obtained from the LGH 
expansion. Since the evaluation of the latter is numeri- 
cally significantly less demanding, we can sample a much 
larger configuration space in this case. To this end, we 
directly enumerate all combinatorially possible ordered 
structures in surface unit-cells of any symmetry and with 
a surface area smaller or equal to a (4 x 4) cell and with 



O coverages up to 1ML. The corresponding data points 
are also shown in Fig. 2J If we first focus on the cover- 
age range up to 0.5 ML, we find the obtained LGH data 
to be fully consistent with the DFT ground state line. 
Namely, there is no structure predicted by the LGH ex- 
pansion that would have an energy that is lower than 
the DFT convex hull, and the LGH Hamiltonian yields 
therefore exactly the same ordered ground states as the 
DFT input data. 



In fact, this is not a coincidental result demonstrating 
the reliability of the achieved LGH expansion, but the 
end product of an iterative procedure. We used the con- 
sistency with the DFT ground state line as another crite- 
rion to judge whether more DFT ordered configurations 
are required as input to the LGH expansion. Initially 
we had started with a smaller number of DFT configura- 
tions as the set discussed above. Having gone through the 
same CV score evaluation, we had identified an optimum 
LGH expansion, but had then obtained LGH data points 
in the direct enumeration leaking below the DFT ground 
state line. Interpreting the corresponding structures as 
important input to the LGH expansion, we would ide- 
ally calculate them with DFT and add them to the DFT 
data set used for the LGH expansion. This was, unfor- 
tunately, not always possible, when the structures pre- 
dicted by the LGH had surface unit-cells that exceeded 
our computational capabilities. In such cases, we looked 
for other structures in smaller unit-cells, which still con- 
tained what we believed were the relevant motifs and 
computed those with DFT. This procedure was repeated 
several times, each time adding new structures to the 
DFT data base, until the present consistent result was 
obtained. 



In the coverage range above ~ 0.5 ML, the situation is 
not that perfect. As apparent from Fig. 0] there are still 
several LGH structures slightly below the DFT ground 
state line. Unfortunately, further improvement along the 
sketched lines is inhibited by the above described propen- 
sity of structures in this coverage range to directly relax 
into geometries with O incorporated below the first Pd 
layer. This renders it very tough to provide new on- 
surfacc O/Pd(100) structures to the data base and im- 
prove on the present LGH expansion. Although not com- 
pletely satisfying, we therefore contend ourselves with the 
achieved expansion. Particular care should therefore be 
exerted, when aiming to use it to describe the higher 
coverage regime, since denser adatom arrangements can 
presumably not be fully reliably described. However, due 
to the overall strongly repulsive interactions, the local oc- 
currence of such denser arrangements at lower coverages 
is rather unlikely in the MC simulations. Correspond- 
ingly, we do expect the results obtained from our expan- 
sion to be reliable for the coverage range below ~ 0.5 ML, 
on which we focus in the present work. 
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FIG. 5: O — T diagram of the critical temperatures T c for the 
order-disorder transition. Shown with stars are the experi- 
mental data from Ref. [Til , and with solid circles the values 
obtained when using the optimum LGH expansion. Addi- 
tionally shown by empty circles are the values obtained when 
using the LDA as exchange-correlation functional in the DFT 
calculations (see text). 



C. Order-disorder transition 

Having established the ground state ordered structures 
we proceed to study the ordering behavior at finite tem- 
peratures. Experimentally, this was investigated in the 
coverage range up to 0.6 ML by Chang and Thiel^. For 
defined initial coverages at the surface, they identified the 
presence of ordered phases at the surface by monitoring 
LEED superstructure spots corresponding to the differ- 
ent periodicities, and the critical temperatures T c (0) for 
the order-disorder transition were determined by the in- 
flection point of the vanishing spot intensities at increas- 
ing temperatures. Avoiding the O-induced reconstruc- 
tion at higher coverages, we focus here on the data in 
the coverage range < 0.35 ML, in which the p(2 x 2) 
or a coexistence of p(2 x 2) and c(2 x 2) phases form 
the ordered structures at low temperatures, cf. Fig. |U 
For this coverage range, Chang and Thiel determined 
the onset of desorption in their ultra-high vacuum ex- 
periments at much higher temperatures than the order- 
disorder transition^ From this we assume that in the 
experiments, the coverage at the surface remained essen- 
tially constant at the initially prepared coverage value for 
all temperatures up to the critical temperatures for the 
order-disorder transition. 

The experimental conditions are simulated by canoni- 
cal MC runs for fixed coverages and at various temper- 
atures. With the definitions in eqs. and ©, our 
order parameters are equivalent to LEED spot intensi- 
ties, so that the determined critical temperatures for the 
order-disorder transition can be directly compared to the 
experimental values. Figure [5] shows the T C (Q) curve ob- 
tained with the optimum LGH expansion together with 
a reproduction of the experimental data. Overall, we 
observe very good agreement, both with respect to the 



absolute temperature values and the trend of increasing 
critical temperatures with coverage. The largest devia- 
tion of about 250 K results at = 0.25 ML, where theory 
predicts a small peak in the critical temperature, which 
is absent in the experimental data and which we dis- 
cuss in Section IVB below. Apart from this feature, the 
agreement with the experimental T C (Q) data is quite sat- 
isfying, if not quantitative. 



D. Population of bridge sites 

The LGH+MC simulations up to now have exclusively 
focused on O adsorption into the fourfold on-surface hol- 
low sites. The already good agreement obtained with ex- 
isting experimental data, together with the significantly 
lower stability of the energetically next favored high- 
symmetry bridge sites apparent in Table HI seem to sug- 
gest that the on-surface low coverage ordering can in- 
deed be understood in terms of the most stable hollow 
sites only. To verify the implied negligible population of 
(possibly metastable) bridge sites, even up to the critical 
temperatures of the order-disorder transition, we proceed 
by including these sites into the LGH expansion. Since 
the intention is at this point only to check on the influ- 
ence of a population of these sites, we consider a reduced 
pool of possible lateral interactions between O atoms ad- 
sorbed in bridge sites, consisting of the equivalents to the 
hollow-hollow pair and trio interactions shown in Fig. [TJ 
Due to the twofold symmetry of the bridge sites, two 
different forms at the same interatomic distance exist for 
some of the interactions, and in addition there is a lateral 
interaction Vq at the very short distance of a/2 between 
O atoms sitting in immediately adjacent bridge sites co- 
ordinated to the same Pd atom. 

However, when computing with DFT configurations 
containing such closely neighboring O atoms at a site 
distance of a/2, we always found them to be unstable 
against relaxation. During the geometry optimization 
the O atoms moved to sites further apart, indicating 
a sizable repulsive Vq. Focusing therefore on 18 or- 
dered configurations that do not contain O adatoms at 
such close distance, the best sets with a varying total 
number of lateral interactions are determined via LOO- 
CV in the same manner as for the hollow-hollow inter- 
actions. Similar to the results in Table HH the differ- 
ent expansions yield consistently the same two dominant 
lateral interactions, namely the first-nearest neighbor 
pair ^(bridge — bridge) and second- nearest neighbor 
pair V£ (bridge — bridge) interactions. Both are largely 
repulsive, with Vy (bridge — bridge) w — 400 me V and 
V 2 P (bridge — bridge) w — 120meV. Turning to the lead- 
ing lateral interactions between O atoms adsorbed in hol- 
low and bridge sites, we again found that structures with 
O atoms in directly adjacent bridge and hollow sites at 
the very short distance of a/y/8 are not stable against 
structural relaxation. As for the other interactions, this 
time the first-nearest neighbor pair interaction and one 
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compact trio interaction turn out to be dominant in the 
LOO-CV based LGH expansion procedure. They are 
also largely repulsive, with values of « — 240 me V and 
—280 meV, respectively. 

We therefore approximately consider the interactions 
involving bridge site O atoms in the MC simulations by 
excluding configurations with O atoms at the above de- 
scribed shortest a/2 and a/v8 distances between bridge- 
bridge and bridge-hollow site pairs, respectively. In ad- 
dition, the two next dominant lateral interactions for 
bridge-bridge and bridge-hollow pairs are explicitly ac- 
counted for using the values stated above. The consec- 
utive MC simulations show virtually no change in the 
ground state ordered structures and critical temperatures 
in the coverage range up to 6 = 0.35 ML. The overall 
largely repulsive interactions together with the signifi- 
cantly less stable on-site energy compared to adsorption 
in the hollow sites, efficiently prevents any significant 
population of bridge sites. For all temperatures up to 
the order-disorder transition, we find less than 10 % of 
the available bridge sites occupied with O atoms, with 
the highest populations obtained for the larger coverages 
in the considered range. To make sure that these results 
are not affected by the uncertainty in the approximately 
determined lateral interactions, we varied the value of 
each of the four lateral interactions by ±100 meV and 
each time reran the MC simulations. This had no influ- 
ence on the findings, so that we do not expect them to 
be invalidated by the crude way of how the bridge-bridge 
and bridge-hollow interactions are considered. Instead, 
we conclude that up to coverages of 8 w 0.35 ML and up 
to the critical temperatures for the order-disorder tran- 
sition, a population of (possibly metastable) bridge sites 
plays no role for the on-surface ordering behavior. 



IV. ACCURACY OF FIRST-PRINCIPLES 
LATERAL INTERACTIONS 



The agreement with the experimental low coverage 
phase diagram (ground state structures and critical tem- 
peratures) suggests that the determined set of first- 
principles lateral interactions is quite reliable. In order 
to get a better understanding of the explicit uncertain- 
ties of the different parameters in this set, we return to 
a critical discussion of all approximations entering the 
LGH approach, and scrutinize their influence on the lat- 
eral interaction values. Uncertainties arise on the one 
hand side due to the truncated LGH expansion and the 
finite number of configurations employed to parameter- 
ize it, and on the other hand due to the approximate 
first-principles energetics, both with respect to total and 
vibrational free energy contributions. 



A. Uncertainties in the LGH expansion procedure 

Table [TT] provides detailed information about the in- 
fluence of most approximations in the LGH expansion 
procedure. Inspecting the basically indistinguishable CV 
score for the expansions with to =9, 10, and 11, one 
might take the scatter in the correspondingly extracted 
lateral interactions as a rough measure for the uncer- 
tainty introduced by truncating the LGH expansion after 
a finite number of terms. Concerning the finite number 
of DFT configurations employed in the parameterization, 
the achieved consistency of the DFT and LGH ground 
state line illustrated in Fig. Q] gives some indication as to 
the minimum number of configurations that is required. 
Correspondingly, the differences in the lateral interaction 
values determined when extending this minimum set by 
two further configurations (upper vs. lower line for each 
expansion in Table |TT| may be taken as reflecting the 
uncertainty due to employing a finite number of DFT 
configurations in the parameterization. 

This leaves as remaining ad hoc feature of our expan- 
sion procedure the decision to not include the on-site 
energy into the fit, but instead fix it at the value of the 
sparsest DFT configuration computed, i.e. the (3 x 3)-0 
structure at 1/9 ML. To this end, we redid all LGH ex- 
pansions in Table [TTI using the same LOO-CV procedure 
described above, but this time including the (3 x 3)-0 
structure into the set of DFT configurations and includ- 
ing the on-site energy £;° n_slte m t the fit. The results 
are remarkably consistent, in the sense that the obtained 
lateral interactions differ in all cases by less than 15meV 
from the values compiled in Table [TTJ and for the expan- 
sions with m = 9, 10, 11 the on-site energy is in fact 
determined at values that are within 15 meV of the com- 
puted binding energy of the (3 x 3)-0 structure. For ex- 
pansions using fewer lateral interaction figures (to < 9) 
this becomes progressively worse, and the increasing in- 
flexibility of the few-interaction expansions starts to shift 
errors between the on-site energy and the lateral inter- 
actions in an uncontrolled way. We therefore conclude 
that for expansions with a sufficient number of interac- 
tion figures it apparently makes very little difference to 
include or not include the on-site energy into the fit; the 
expansions are stable in this respect. In view of the sig- 
nificantly different inaccuracies in the determination of 
the on-site energy and lateral interactions discussed be- 
low, we nevertheless prefer to fix the on-site energy at 
the value of the sparsest DFT configuration computed. 

Summarizing the discussion on the uncertainties in the 
LGH expansion and parameterization procedure, we es- 
timate the uncertainty introduced by the various approx- 
imations to be of the order of 10-20 meV in the dominant 
lateral interactions. When using the first- principles lat- 
eral interactions to determine quantities describing the 
mesoscopic ordering behavior, this uncertainty needs to 
be taken into account. Figure [5] illustrates this for the 
determined critical temperatures for the order-disorder 
transition by comparing the results obtained for the dif- 
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eq. (U) solely to the NoEf, of the different computed 
configurations. The formally correct procedure would be 
to evaluate for each configuration also the vibrational 
free energy contribution to the binding energy, i^ lb " (T) , 
and consider it in the parameterization. As apparent 
from eq. ©, what determines this vibrational contri- 
bution are the changes of the vibrational modes during 
adsorption. Since this will predominantly concern the 
adsorbate-derived vibrational modes, we estimate the or- 
der of magnitude of this contribution by the zero point 
energy (E ZPE ) difference arising from the change of the 
~02 ' 015 ® 2 s ^ re ^ cn frequency, wq 2 ( ga s) > to the O-substrate stretch 



FIG. 6: O — T diagram of the critical temperatures T c for the 
order-disorder transition. Compared are the curves obtained 
when using the first-principles lateral interactions for different 
LGH expansions compiled in Table ITT1 and each time using 24 
DFT configurations in the parameterization: Solid line, full 
circles: m = 9; dashed line, empty circles: m = 10; dotted 
line, empty circles: m = 11. 



ferent LGH expansions with m =9, 10, and 11 in Table 
im for which the extracted lateral interactions exhibit a 
scatter of about the estimated magnitude. The critical 
temperatures show a variation of less than 100 K and the 
overall form of the T c (0)-curve is almost untouched. The 
systematics of the present LGH expansion procedure pro- 
vides thus an error-controlled link between the electronic 
structure calculations and the mesoscopic statistical sim- 
ulations, which allows to assess which quantities can be 
determined with which uncertainty for a given accuracy 
of the underlying first-principles energetics. 



B. Uncertainties in the first-principles energetics 

The energetic parameters in the LGH expansion in eq. 
(fTj) comprise total and vibrational free energy contribu- 
tions. Uncertainties in our approach arise therefore out of 
the treatment of the vibrational free energy contribution 
and the approximate DFT energetics, where the latter 
contains uncertainties due to the numerical setup and 
due to the approximate exchange-correlation functional. 
We will discuss these three sources of uncertainties sub- 
sequently. 



mode, i.e. we approximate eq. ([6]) roughly with 
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where u>q /p<j(ioo) (*) is the O-substrate stretch frequency 
of the ith adsorbed O atom in the corresponding con- 
figuration. Within the frozen substrate approxima- 
tion we calculate the latter stretch mode in good 
agreement to the experimental data^ as hu) p (2x2) = 
50meV in the p(2 x 2)-0 configuration and hu) c (2x2) = 
33meV in the c(2 x 2)-0 configuration. Compared 
to the computed stretch frequency of ?ia;o 2 (gas) = 
190 meV of the O2 molecule, this yields an esti- 
mated vibrational contribution to the binding energy of 
-(^ P (2x2)-^o 2 (gas)/2)/2 = 23meV and -(?kJ c (2x2)- 
?kJo 2 (gas)/2)/2 = 31meV per O atom for the two config- 
urations, respectively. In the LGH expansion these con- 
tributions get separated into a non-coverage dependent 
part, which enters the on-site energy, and a coverage- 
dependent part, which enters the lateral interactions. 
Taking the difference between the estimated p(2 x 2)-0 
and the c(2 x 2)-0 vibrational free energy contributions 
as a measure of the coverage dependence we thus arrive at 
an uncertainty of ~ 5 meV in the lateral interactions and 
an uncertainty of ~ 30 meV in the on-site energy due to 
the neglect of vibrational contributions in our approach. 
The on-site energy is thus much more affected by the pre- 
dominantly non-coverage dependent shift in vibrational 
frequency from the O2 stretch to the O-substrate stretch 
mode. 



Numerical uncertainties in the DFT energetics 



1. Vibrational free energy contribution 

Separating each of the energetic parameters in the 
LGH (i.e F° n ~ sitc , VP, V£, ...) into a total energy term 
and a term due to the vibrational free energy contribu- 
tion, one arrives at a LGH expansion for the total bind- 
ing energies and a LGH expansion for the vibrational free 
binding energies. In Section IIIB the vibrational part was 
completely neglected by equating the right hand side of 



Turning to the effect of the approximate DFT total 
energies, we first distinguish between numerical inac- 
curacies which arise out of the finite basis set, the k- 
point sampling or the chosen supercell geometries, and 
the more fundamental uncertainty due to the employed 
approximate xc functional. In principle, the numerical 
inaccuracies can be reduced to whatever limit desired, 
although this may quickly lead to unfeasible computa- 
tions in practice. In this respect, the finite slab and vac- 



11 



uum thicknesses employed in our supercell setup are the 
end result of extensive test calculations, in which these 
quantities were progressively increased until the relative 
binding energies were converged to the desired limit of 
±10meV/O atom. An important point to note concern- 
ing the k-point sampling is that compatible k-meshes 
must be employed in calculations involving different sur- 
face unit-cells. Compatible in this respect means that al- 
ways the same points within the Brillouin zone are sam- 
pled, even when the latter changes its size due to the 
different periodicity of the real-space surface unit-cell. 
As long as this is considered, we found only a negligi- 
ble influence on the LGH parameters, when increasing 
the k-point density to higher values than the one of our 
standard computational setup described in Section IIC. 
This leaves as most notable numerical approximation the 
finite energy cutoff for the plane wave representation in 
the interstitial region between the muffin tin spheres in 
the LAPW/APW+lo scheme. To this end, we repeated 
the LGH expansion with the optimum m — 9 lateral in- 
teraction figure set, cf. Table [H] using as input the 24 
DFT configurations with total energies computed self- 
consistently at cutoffs of 18 Ry, 20 Ry and 22 Ry (20 Ry 
corresponding to our standard basis set size) . The varia- 
tions in the determined values for all lateral pair, trio and 
quattro interactions are within 5meV. In contrast, what 
changes notably is the on-site energy (determined by the 
(3 x 3)-0 structure). With increasing cutoff, its value 
decreases by more than 100 meV, and from the conver- 
gence pattern using energy cutoffs up to 30 Ry we expect 
it to decrease by another ~ 50 meV, when extrapolating 
to infinite cutoff. The reason for this different conver- 
gence behavior is again that the on-site energy gathers 
all non-coverage dependent contributions. As such, it in- 
cludes also all the inaccuracies in the description of the 
isolated O2 molecule, the total energy of which enters 
the binding energy of all DFT configurations used in the 
parameterization alike. The slow convergence of the on- 
site energy is therefore mainly caused by the known slow 
convergence in the description of the gas phase molecule. 
At feasible basis set sizes (for full potential calculations 
of the present system), the lateral interactions can there- 
fore again be determined at a much higher accuracy than 
the on-site energy. 



3. Approximate exchange-correlation functional 

The uncertainty introduced by the approximate xc 
functional can not be determined in a comparably quanti- 
tative manner as for the numerical approximations, since 
the exact functional is not known. In order to get at least 
a feeling for the scatter in the results when using differ- 
ent present-day xc functionals, we evaluated the lateral 
interactions also within the local density approximation 
(LDA) 32 , which is a functional that by its very construc- 
tion is known to yield significantly different adsorbate 
binding energies than the hitherto employed GGA-PBE 



functional. Most prominently, the two functionals yield 
e.g. binding energies for the free O2 molecule that differ 
by 1.35 eV (!) when computed with our accurate full- 
potential LAPW/APW+lo scheme. All 24 DFT input 
configurations were correspondingly calculated fully re- 
laxed using the LDA and the LDA optimized lattice con- 
stant (a = 3.836 A). Using their energetics in the LGH 
expansion procedure, we obtain the lateral interaction en- 
ergies compiled in Table IIIII for the previously discussed 
optimum set with m = 9 interaction figures. 

Comparing to the results obtained before using the 
GGA-PBE as xc functional, also compiled in Table IIII| 
one clearly observes a similar pattern as before for the 
numerical uncertainties in the sense that the difference 
between the two xc functionals shows up predominantly 
in the on-site energy. In fact, the lateral pair, trio and 
quattro interactions are strikingly little affected, consid- 
ering that the two xc functionals yield adsorbate binding 
energies that differ on the order of ~ 0.7 eV. The reason 
is again as before, i.e. most of this variation arises out 
of the description of the free gas phase molecule, which 
affects all computed configurations alike and therefore 
solely enters the non-coverage dependent LGH parameter 
^.on-site ^ rpj^ j cnown i ar g er variation of adsorbate bind- 
ing energies when using different present-day xc function- 
als is often put forward as a generic argument against the 
reliability or usefulness of first-principles lateral interac- 
tions. At least for the present system, we can now qualify 
this concern. Our analysis shows that the on-site energy 
can indeed not be determined with a high accuracy, cf. 
Table IIIII However, the lateral interactions themselves 
can be. For the O/Pd(100) system, the uncertainty due 
to the approximate xc functional seems to be of the order 
of 50meV for the dominant pair interactions. The far- 
ther ranging and higher-order interaction terms in Table 
IIIII exhibit intriguingly an even smaller scatter and we 
are currently pursuing further studies to elucidate the 
generality of this result. 

The largely different accuracy with which on-site en- 
ergy and lateral interactions can be determined allows 
also to scrutinize when the first-principles parameters 
may be employed to reliably describe mesoscopic system 
quantities. We focus here on the effect of the approxi- 
mate xc functional, since this leads to a larger uncertainty 
compared to the discussed treatment of the vibrational 
free energy contributions and the numerical inaccuracies 
in the total energies. Since the on-site energy has no ef- 
fect on the order-disorder transition in a canonic adsor- 
bate ensemble, the small variation in the lateral interac- 
tion energies exhibited in Table IIIII rationalizes the good 
agreement with the experimental critical temperatures. 
Using the determined LDA interaction values we indeed 
obtain a T C (Q) curve that is very similar to GGA-PBE as 
shown in Fig. [5j The overall shape of the curve includ- 
ing the peak structure is completely preserved, and the 
absolute T c values are merely shifted by up to 200 K to 
higher temperatures. This variation may thus be taken 
as an indication of the accuracy with which this quan- 
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TABLE III: First-principles lateral interactions obtained using the LDA and the GGA-PBE as exchange-correlation functional 
for the optimum set of m = 9 interaction figures determined in Table [II] Additionally shown are the values obtained, when 
restricting the LGH expansion in eq. ((!]) to pairwise interactions only. Units are meV. 
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tity can presently be calculated from first-principles. For 
other quantities, where the on-site energy enters explic- 
itly, this accuracy will be significantly worse. Staying 
with the presently available experimental data for the 
O/Pd(100) system^, this would e.g. be the case when 
aiming to compute the measured temperatures for oxy- 
gen desorption, where the on-site energy enters explicitly. 

One can, of course, always argue in terms of a smaller 
uncertainty by favoring a particular approximate xc func- 
tional over others (e.g. because the functional is known 
to reproduce well similar system quantities). In the sense 
of first-principles theory, the more appropriate approach 
would, however, be to eliminate errors by using higher- 
level theory. In this respect, one could suitably exploit a 
finding like for the present system that most of the xc un- 
certainty gets concentrated in the on-site energy by e.g. 
computing only the latter with a higher-level technique 
(or appropriate xc-correction scheme^) and using the al- 
ready quite accurate lateral interactions determined by 
present-day xc functionals. 

V. COMPARISON TO EMPIRICAL 
PARAMETERS 

Instead of a first-principles determination, the tradi- 
tional method to obtain the strengths of lateral inter- 
actions has been to adjust the predictions from atom- 
istic lattice-gas models to experimental observations of 
the coverage and temperature dependence of on-surface 
ordering. To keep the number of free fit parameters 
low, the focus has often been on pair interactions only, 
i.e. eq. fl} was truncated after the pair terms, and 
usually also only the dominant short-ranged pair in- 
teractions were considered. In this respect, extensive 
model studies££2&±&±ii±^ have shown that the p(2 x 2) 
and c(2 x 2) ordering frequently found at (100) cubic 
metal surfaces requires nearest neighbor repulsions that 
are so strong that they effectively yield a site exclusion 
(—Vi >> ksT), as well as weaker second nearest neigh- 
bor repulsions. In order to produce a true p(2 x 2) order- 
ing further interactions are required. This can either be 
a third nearest neighbor attraction, or a fourth nearest 
neighbor repulsion, or more generally a combination of 
these interactions that fulfills V£ -V£/2> 

Comparing with our first-principles data compiled in 



Table [TTTI we see that they perfectly fit into this expected 
qualitative pattern: Strong first nearest neighbor repul- 
sion, weaker second nearest neighbor repulsion, and the 
combination V£ — Vf /2 is attractive. However, differ- 
ences arise when turning to a more detailed modeling 
of the experimental O/Pd(100) phase diagram of Chang 
and Thieli 4 . For a lattice-gas model considering nearest- 
neighbor site exclusion and finite second and third near- 
est neighbor interactions, Liu and Evans showed that the 
best overall topological agreement with the experimental 
phase diagram is achieved for V% ~ — VJ?*^ In particular, 
the position of the disordered to c(2 x 2)-0 transition line 
shifts notably with the ratio of —V^/Vq, and for a ratio 
of ~ 1 the T C (Q) curve no longer exhibits the peak struc- 
ture at O = 0.25 ML visible in the first-principles curve 
in Fig. [5] Since this resembles the shape of the experi- 
mental curve better, it could indicate that the calculated 
first-principles ratio of —V£/V£ ~ 3 is too large. On the 
other hand, it is intriguing to see that both LDA and 
GGA-PBE yield roughly the same —V^/V^ ratio and in 
turn clearly a peak in the T c (0) curve, cf. Fig. [5] Since 
Chang and Thiel discuss also a rather large uncertainty of 
~ 0.05 ML in their experimental coverage calibration^, 
a careful remeasurement of the experimental T C (Q) curve 
would thus be very interesting to fully settle this point. 

A much more fundamental difference between the set of 
empirical and first-principles lateral interactions is that 
the former is restricted to just pairwise interactions. To 
test the validity of this frequently applied approxima- 
tion, we restrict our LGH expansion to pair interactions 
up to third nearest neighbors and repeat the parameter- 
ization using exactly the same 24 DFT configurations as 
in Section IIIB. The resulting interactions are also com- 
piled in Table |In| using the LDA or GGA-PBE as xc 
functional. Comparing with the values obtained for each 
functional from the proper LGH expansion (and using 
exactly the same first-principles input), we note signifi- 
cant differences. The dominant short-ranged interactions 
are overestimated, in the case of the first nearest neigh- 
bor pair interaction by ~ 60mcV. Overall, the lateral 
interactions are too repulsive, which is a consequence of 
them having to mimic the effect of the neglected overall 
repulsive trio interactions. As expected, the true micro- 
scopic lateral interactions are thus blurred by spurious 
contributions arising from the invalid assumption that 
lateral interactions at metal surfaces could be expressible 
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FIG. 7: O — T diagram of the critical temperatures T c for the 
order-disorder transition. Shown with stars are the experi- 
mental data from Ref. [Til , and with solid circles the values 
obtained when using the optimum GGA-based LGH expan- 
sion, cf. Table Mil These curves are compared to those ob- 
tained, when considering only pair interactions in the LGH: 
Using the GGA values compiled in Table Mil (empty circles) , 
or the empirically adjusted values by Liu and Evans (empty 
squares), see text. 



0.1 ML outside the fitted coverage range, the predicted 
critical temperature for the order-disorder transition is 
wrong by 1200 K. 

The spurious contributions contained in the effective 
pair interactions also prohibit a systematic determina- 
tion of which pairwise interactions to consider. In the 
systematic LGH expansion detailed in Section IIIB, the 
lateral interactions are chosen out of a larger pool of lat- 
eral interactions and the ultimately determined pair in- 
teractions decrease in magnitude with distance. In con- 
trast, when we extend the pair-only expansion to include 
the fourth and fifth nearest neighbor interaction, we ob- 
tain a negligible of -2meV, but a V 5 P of -18meV. 
The spurious contributions contained in the effective V 5 P 
thus feign that such far ranging interactions are signifi- 
cant, whereas the proper first-principles LGH expansion 
shows no significant interactions between adsorbates at 
distances larger than the fourth nearest neighbor site. In 
empirical approaches this uncertainty (or degree of ar- 
bitrariness) as to which interactions to consider is even 
further aggravated, since there are typically several sets 
of pairwise interactions which equally well fit the less 
stringent mesoscopic ordering behavior. 



in pair-only terms. The invalidity of this assumption is 
already clearly demonstrated by the similar magnitude of 
the dominant trio interactions (V*, Vj, etc. in Tabic |LH|) 
compared to the dominant short-ranged pair interactions, 
which necessarily reduces pair-only restricted expansions 
to a set of effective parameters at best. 

The deficiencies introduced by such an effective de- 
scription are not only the lack of microscopic meaning of 
the lateral interactions themselves, but also errors in the 
mesoscopic system properties calculated with these inter- 
actions. This is again exemplified for the order-disorder 
transition in Fig. [JJ which shows that the T C (Q) curve 
derived from the pair-only lateral interactions in Table 
IIIII differs notably from the curve derived from the corre- 
sponding proper LGH expanded lateral interactions. The 
critical temperatures are up to 350 K higher, even though 
being based on exactly the same first-principles energet- 
ics. The shortcomings due to the restricted LGH can also 
not be overcome, when turning to a completely empiri- 
cal approach and adjusting the lateral interaction values 
to reproduce experiment. Within the above described 
lattice-gas model with pair interactions to third nearest 
neighbor and fitting to the critical temperatures in the 
coverage range below 0.25 ML, Liu and Evans determined 
— — V 2 P w 40meVA2., i.e. lateral interactions that dif- 
fer considerably from the first-principles values compiled 
in Table IIIII As apparent from Fig. [7j these empirical 
parameters do indeed lead to an excellent fit of the exper- 
imental T c {&) curve in the coverage range up to 0.25 ML. 
However, outside the fitted coverage range, the effective 
description rapidly breaks down and predicts a critical 
temperature curve in qualitative disagreement with ex- 
periment. Already at a coverage of 0.35 ML, i.e. only 



VI. CONCLUSIONS 

In summary, we have presented a first-principles 
lattice-gas Hamiltonian study of the on-surface ordering 
behavior of O adatoms at the Pd(100) surface. A key 
feature of the approach is the systematics of the LGH 
expansion, both with respect to the selection of the con- 
sidered lateral interactions and with respect to the or- 
dered configurations, the DFT energetics of which is em- 
ployed in the parameterization. In contrast to empirical 
or semi-empirical ad hoc parameterizations the approach 
does provide an error-controlled link between the elec- 
tronic structure regime and the system's mesoscopic en- 
semble properties. Carefully scrutinizing the approxima- 
tions entering at the level of the underlying DFT ener- 
getics, at the LGH interface and in the statistical simula- 
tions, we identify the approximate exchange-correlation 
functional employed to obtain the DFT energetics as the 
major source of uncertainty in the approach. This uncer- 
tainty affects predominantly the on-site energy, and only 
to a much lesser degree the lateral interactions them- 
selves. Comparing LDA and GGA exchange-correlation 
functionals, we estimate the accuracy of the latter to be 
within ~ 60 meV for the present system. This shows 
that the known much larger variation of adsorbate bind- 
ing energies with different xc functionals can not be used 
to uncritically question the reliability or usefulness of lat- 
eral interactions derived from first-principles per se. 

The rather high accuracy of the first-principles lateral 
interactions also rationalizes the obtained good agree- 
ment with the experimental low coverage phase diagram 
for the O/Pd(100) system, i.e. the ground state or- 
dered structures and critical temperatures for the order- 
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disorder transition. An important feature of the set of 
lateral interactions is that it contains many-body inter- 
actions of comparable magnitude to the dominant short- 
ranged pair interactions. This demonstrates the invalid- 
ity of the assumption of exclusively pairwise interactions 
between adsorbates at metal surfaces that is frequently 
applied in empirical approaches. Indeed, when restricting 
the LGH expansion to just pairwise lateral interactions, 
we find that this results in effective interactions which 
contain spurious contributions that are of equal size, if 
not larger than any of the uncertainties e.g. due to the 
approximate xc functional. These effective parameters 
lack microscopic meaning and lead to uncontrolled er- 
rors in the mesoscopic system properties calculated with 
them. In the present system, this is illustrated clearly by 
critical temperatures that are up to 350 K higher than 
the ones obtained with the proper set of lateral interac- 
tions, even though both sets were based on exactly the 



same first-principles energetics. 

We conclude that at least for the present system, an 
appropriate first-principles statistical mechanics frame- 
work and present-day DFT energetics can determine lat- 
eral interactions at least as reliably as the traditional em- 
pirical approaches, and without suffering from the non- 
uniqueness in the selection of pairwise-additive adspecics 
interactions which reasonably fit available data. 
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